function diff = KFE_t(m_old,Q_t,delta,param,aux)

    % Create new and old grids for m
    lnm=ln_m_KFE([param.m_l,param.m_u],1,aux);
    lnm_old = ln_m_KFE(m_old,1,aux);

    % Assign drift taking into account the difference in grid points
    mu=param.mu-(lnm-lnm_old)./aux.dt;

    % Calculate derivative of delta
    delta_p=[[delta(2:end)-delta(1:end-1)]./[lnm(2:end)-lnm(1:end-1)]; nan(1)];


    %% Create useful Matrices

    % Level (diagonal with ones)
    T_l=speye(length(lnm));

    % First derivative (first derivative, central, backward and forward)
    [T_p, T_l_p, T_u_p]=M_p(lnm,ones(size(lnm)));

    % Second derivative
    T_pp=M_pp(lnm);

    % Separation threshold q_l = q(m_l)
    T_q_l=sparse(aux.n_m,aux.n_m);
    T_q_l(:,1)=1;

    % Derivative at separation theshold q_p_l = q'(m_l)=
    T_q_p_l=sparse(aux.n_m,aux.n_m);
    T_q_p_l(:,1)=-1./(lnm(2)-lnm(1));
    T_q_p_l(:,2)=1./(lnm(2)-lnm(1));

    if param.K>0

        % Indicator for replacement region q(m_h) 
        ind_r=find(exp(lnm)>param.m_h,1)-1;
        T_q_r=sparse(aux.n_m,aux.n_m);
        T_q_r(:,ind_r)=1;

        % Central derivative at replacement region q_p_r = q\prime(m_h)
        T_q_p_r=sparse(aux.n_m,aux.n_m);
        T_q_p_r(:,ind_r-1:ind_r+1)=repmat(T_p(ind_r,ind_r-1:ind_r+1),aux.n_m,1);

        % Second derivative at replacement region T_q_pp_r = q\prime\prime(m_h)
        T_q_pp_r=sparse(aux.n_m,aux.n_m);
        T_q_pp_r(:,ind_r-1:ind_r+1)=repmat(T_pp(ind_r,ind_r-1:ind_r+1),aux.n_m,1);
    end

    % Indicator for the tio of the support q_p_u =  q(m_u)
    T_q_u=sparse(aux.n_m,aux.n_m);
    T_q_u(:,end)=1;

    %% Create transition matricies

    if param.K>0
        T=(param.zeta+delta+(1-param.alpha).*delta_p).*T_l ...
            +(mu-param.sigma^2/2+(1-param.alpha).*(param.zeta+delta)).*T_p ...
            -param.sigma^2/2.*T_pp...
            -(1-1./param.s).*(param.zeta+delta(1)).*T_q_l ...
            +(1-1./param.s).*param.sigma^2/2.*1./(1-param.alpha).*T_q_p_l ...
            -1./param.s.*param.zeta.*T_q_u ;
    else 
        T=(param.zeta+delta+(1-param.alpha).*delta_p).*T_l ...
            +(mu-param.sigma^2/2+(1-param.alpha).*(param.zeta+delta)).*T_p ...
            -param.sigma^2/2.*T_pp...
            -(1-1./param.s).*(param.zeta+delta(1)).*T_q_l ...
            +(1-1./param.s).*param.sigma^2/2.*1./(1-param.alpha).*T_q_p_l ...
            -1./param.s.*param.zeta.*T_q_u ;    
    end

    % Impose beverage curve (inflow and outflow from unemployment)
    T(1,:)=-1./param.s.*(param.zeta.*T_q_u(1,:)     +param.sigma^2./2./(1-param.alpha).*T_q_p_l(1,:)-(param.zeta+delta(1)).*T_q_l(1,:)); 

    % Adjust transition matrix for timeperiod dt 
    A=speye(aux.n_m)+aux.dt.*T;
    % Impose that the beverage curve holds at the top
    A(end,:)=[zeros(1,aux.n_m-1),1]+[A(1,:)-[1,zeros(1,aux.n_m-1)]].*(1-param.s);
    % Use implicit method
    Ainv=A^(-1);

    % Impose "smooth pasting" at the m_l
    a=((1-param.alpha).*(param.zeta+delta(1)))./(lnm(2)-lnm(1));
    b=param.sigma^2./2.*1./(lnm(2)-lnm(1))./((lnm(3)-lnm(1))/2);
    c=param.sigma^2./2.*1./(lnm(3)-lnm(2))./((lnm(3)-lnm(1))/2);
    Ainv(2,:)=c./(a+b+c).*Ainv(3,:)+(a+b)./(a+b+c).*Ainv(1,:);

    diff=Ainv^(aux.n)*Q_t;

end